function [X, L, X_err, L_err, iter, flag] = invit(A, Xo, X_ref, L_ref, maxit, tol, print)
%
%  [X,L,X_err,L_err,iter,flag] = invit(A, Xo, X_ref, L_ref, maxit, tol, print)
%
%  INVIT applies inverse iteration to the eigenvalue problem
%
%    (A - L*I)*X = 0
%
%  The iterations are seed by a handful of power iterations to give a fair
%  estimate of the dominant eigenvalue.  We know this should be about 1, but for
%  early values of keff, there may be other values above 1, so e'd like to avoid
%  catching one of those.
%
%    X[0]   = normalized guess, i.e. ||X[0]||_2 = 1.0
%    Solve (A - mu*I) & X[l+1] = X[l]
%    X[l+1] = X[l+1] / ||X[l+1]||_2
%    L[l+1] = X[l+1]' * A * X[l+1]
%
%  To solve, we use Matlab's built in GMRES without preconditioning.
%
%  Inputs:
%
%    real       A(N,N)     -- N-by-N matrix
%    real       xo(N)      -- Initial guess N-vector.
%    real       lambda_ref -- The actual eigenvalue for checking errors
%    integer    maxit      -- The maximum number of iterations.
%    real       tol        -- The error tolerance (defined...?)
%
%  Outputs:
%
%    real       x(N)       -- The computed eigenvector.
%    real       lambda     -- The computed eigenvalue.
%    real       rho        -- The estimated dominance ratio.
%    real       err_norm   -- The norm of the error.
%    integer    iter       -- The number of iterations performed.
%    integer    flag       -- flag (0=okay, 1=maxit reached, -1=bad bad news)
%

if (nargin==6)
    print = 1;
end
flag    = 1; % unconverged
X_err   = 0;
L_err   = 0;
L1 = 0;
% Precondition the initial guess with a few power iterations
[X L]  = powermethodRHO(A,Xo,10);
scale  = 1;
% iA      = A - L*speye(length(A));
setup.droptol = 0.01;
setup.type    = 'ilutp';
% [aL aU] = ilu(iA);%,setup);
% redo = 1;
L = 1; % know L is close to 1. Make this our first guess.
for iter = 1:maxit
    
    % compute inverse operator (A-mu*I)
    toll = 1/(scale*iter*10);
    disp(sprintf('%6.6e, %6.6e',abs(L-L1),toll));
    
    %   if (abs(L-L1) > toll && iter > 1 && redo == 1)
    if (iter < 2)
        iA         = A - L*speye(length(A));
        setup.droptol = 0.001; % 0.005 - 0.01 works
        disp( 'recomputed ilu ' )
        [aL aU ]   = ilu(iA,setup);
        redo = 0;
    end
    % compute ILU(0) preconditioner
    
    [X flag relres it] = gmres(iA,X,20,1e-7,100,aL,aU,X*scale);
    fprintf(' gmres flag = %4i \n ', flag )
    % save normalization for next time and then scale
    scale = norm(X);
    X     = X / scale;
    
    % update L and err est's
    L1    = L;
    L     = X'*(A*X);
    X_err = norm(X - sign(X(1)/X_ref(1))*X_ref, 2);
    L_err = abs(L - L_ref);
    if ( X_err <= tol  && L_err <= tol)
        flag = 0;
        break;
    end;
    printout(iter,it(1),it(2),X_err,L,L_err,print)
end
if (flag==0 && print == 1)
    disp(' *** InvIt: Final Results *** ');
end
if (flag==1 && print==1)
    disp(' *** Warning: InvIt did not converge ***');
end
printout(iter,it(1),it(2),X_err,L,L_err,print)

end

function printout(iter,i2,i3,X_err,L,L_err,print)
if (print==1)
    disp(sprintf(...
        ' Iter = %5i %5i %5i, X_err = %6.5e, L = %6.4e, L_err = %6.3e', iter, i2,i3, X_err,L,L_err));
end
end

